# -*- coding: utf-8 -*-
"""
Created on Fri Jun 07 17:15:15 2013

@author: Maxim
"""
from airfoil import Airfoil
from propeller import propeller
from numpy import array,linspace,zeros,polyfit,poly1d,vstack,transpose,array_equal
from scipy.optimize import minimize,fsolve
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
import engine

def _interpolate(x,y,xnew,method):
    if method=='spline':
        curve = interp1d(x,y,'cubic')
    elif method=='poly':
        a = polyfit(x,y,len(x)-1)
        curve = poly1d(a)
    ynew = curve(xnew)
    return ynew


def create_propeller(x,afOld):
    diameter = 1.75
    hubDiameter = 0.3
    hubX = hubDiameter / diameter
    baseAfPath = 'clark-Y.txt'

    chordX = array([hubX,0.35,0.5,0.75,0.95])
    pitchX = array([hubX,0.35,0.5,0.75,0.95])
    thickX = array([hubX,0.35,0.5,0.75,0.95])
    xnew = linspace(hubX,0.95,10)
    betaSet = pitchX[3]
    n1,n2,n3 = len(chordX), len(pitchX), len(thickX)
    chordY = x[0:n1]
    pitchY = x[n1:n1+n2]
    thickY = x[n1+n2:]

    method = 'spline'
    chord = _interpolate(chordX,chordY,xnew,method)
    pitch = _interpolate(pitchX,pitchY,xnew,method)
    thick = _interpolate(thickX,thickY,xnew,method)
    
    prop = propeller()
    if afOld==None:
        af = list()
        for tc in thick:
            afnew = Airfoil()
            afnew.read_txt(baseAfPath)
            afnew.analyze_geometry()
            afnew.scale_thickness(tc,False)
            afnew.build_aero_table(alphaSeq=[-40.0,40.0,1.0])
            af.append(afnew)
        prop.airfoil = af
    else:
        prop.airfoil = afOld

    prop.name = 'MDO_project'
    prop.diameter = diameter
    prop.hubDiameter = hubDiameter
    prop.chord = chord
    prop.beta = pitch
    prop.x = xnew
    prop.r = xnew*diameter/2.0
    prop.radius = diameter/2.0
    prop.numBlades = 3.0
    prop.thickness = thick
    prop.numStations = len(xnew)
    prop.analyze_geometry()
    prop.set_beta(betaSet)
    
#    if show:
#        fig = plt.figure(1)
#        ax1 = fig.add_subplot(311)
#        ax1.grid(True)
#        ax1.hold(True)
#        ax1.plot(xnew,chord)
#        ax2 = fig.add_subplot(312)
#        ax2.grid(True)
#        ax2.plot(xnew,pitch)
#        ax3 = fig.add_subplot(313)
#        ax3.grid(True)
#        ax3.plot(xnew,thick)
#        plt.show()
    return prop


class FixedPitchThrust:
    def __init__(self,engine):
        self.nFeval = 0
        self.niter = 0
        self.engine = engine
        self.prop = None
        self.thicknessOld = zeros(3)
        self.S = 11.4
        self.W = 620.0*9.81
        self.CD0 = 0.04
        self.k = 0.03985
    def update_propeller(self,x):
        self.prop = create_propeller(x)
    def _vel_obj(self,x,pset,rho,weight):
        self.nFeval += 1
        rpm, V = x[0],x[1]
        beta = self.beta
        data = self.prop.analyze_prop(beta,rpm,V,rho)
        Tav = data[0]
        Preq = data[1]
        Pav = self.engine.get_availablePower(rho,rpm,pset)
        q = rho*V*V/2.0
        qS = q*self.S
        CL = weight /qS
        CD = self.CD0 + self.k*CL*CL
        Treq = CD*qS
        return [Tav-Treq,Pav-Preq]
    def velocity_at_power_setting(self,pset,rho,weight):
        x = fsolve(self._vel_obj,[1000.0,40.0],args=(pset,rho,weight),xtol=1e-10)
        rpm, V = x[0],x[1]
        return rpm,V
    def max_efficiency(self,pset,rho):
        x0 = array([0.0413,0.0903,0.1135,0.0921,0.0572,
                44.0,37.5,29.5,23.0,20.0,
                0.4,0.24,0.1175,0.0987,0.0822])
        lb = 0.75*x0
        ub = 1.25*x0
        def objective(x,pset,rho):
            thicknessNew = array(x[-5:])
            if array_equal(thicknessNew,self.thicknessOld):
                afOld = self.prop.airfoil
            else:
                afOld = None
            self.thicknessOld = thicknessNew
            self.prop = create_propeller(x,afOld)
            self.beta = x[8]
            #out = 'Efficiency: '
            #SAR_obj = 0.0
            effy = 0.0
            for p in pset:
                for w in self.W:
                    rpm,V = self.velocity_at_power_setting(p,rho,w)
                    data = self.prop.analyze_prop(self.beta,rpm,V,rho)
                    effy += data[-1]
#                    ff = self.engine.get_fuelFlow(rho,rpm,p)
#                    SAR = V/ff/1000.0
#                    out += '\t%.4f'%(SAR)
#                    SAR_obj += SAR
#            out += '| %.4f'%SAR_obj
            #print out
            out = effy / 6.0
            print out
            return -out
        bnd = transpose(vstack([lb,ub]))
        bnds = tuple(tuple(x) for x in bnd)
        args = (pset,rho)
        xOpt = minimize(objective,x0,bounds=bnds,args=args,method='SLSQP')
        print xOpt
        print self.nFeval
        print self.niter
        newProp = create_propeller(xOpt.x,None)
        newProp.plot_propeller()


def run_design():
    eng = engine.load('ROTAX_912ULS')
    thrust = FixedPitchThrust(eng)
    thrust.W = array([550.0,620.0])*9.81
    pset = array([50.0,75.0,95.0])
    thrust.max_efficiency(pset,1.2255)

def save_designs():
    x0 = array([0.0413,0.0903,0.1135,0.0921,0.0572,
                40.0,33.5,25.5,16.0,
                0.4,0.24,0.1175,0.0987,0.0822,19.0])
    xOpt = array([  3.09750000e-02,   6.77250000e-02,   8.51250000e-02,
         6.90750000e-02,   4.29000000e-02,   4.99793255e+01,
         3.67684474e+01,   2.86244787e+01,   1.32814750e+01,
         3.99968568e-01,   2.26012901e-01,   1.14543646e-01,
         9.21425537e-02,   7.59959892e-02,   2.85043362e+01])
    
    prop0 = create_propeller(x0,None)
    propOpt = create_propeller(xOpt,None)
    prop0.betaRange = [15.0,30.0]
    propOpt.betaRange = [15.0,30.0]
    prop0.name = 'MDO_initial'
    propOpt.name = 'MDO_optimum'
    for i in range(len(prop0.airfoil)):
        name = prop0.airfoil[i].name
        prop0.airfoilName.append(name)
    for i in range(len(propOpt.airfoil)):
        name = propOpt.airfoil[i].name
        propOpt.airfoilName.append(name)
    prop0.write_xls()
    propOpt.write_xls()

def get_required_data():
    eng = engine.load('ROTAX_912ULS')
    prop = propeller()
    prop.read_xls('MDO_initial')
    beta = 19.0

    thrust = FixedPitchThrust(eng)
    thrust.beta = beta
    thrust.prop = prop
    W = 620.0*9.81
    rho = 1.2255
    pset = array([50.0,75.0,95.0])
    vel = list()
    LD = list()
    effy = list()
    ff = list()
    sfc = list()
    print 'pset\t fuelFlow kg/s\tV,m/s\teffy\tP,W'
    for p in pset:
        rpm,V = thrust.velocity_at_power_setting(p,rho,W)
        data = prop.analyze_prop(beta,rpm,V,rho)
        effy.append(data[-1])
        vel.append(V)
        fuelFlow = eng.get_fuelFlow(rho,rpm,p)
        ff.append(fuelFlow)
        print '%.1f\t%.2e\t%.4f\t%.4f\t%.4f'%(p,fuelFlow,V,data[-1],data[1])
        

if __name__=="__main__":
    run_design()